######################Simulations for posterior predictive p-values#############
###################################plot function################################

dta_plot <- function(ppp,p.values=NA,ppp_boot=NA,bootppp=NA,m='single'){
  library(tidyverse)
  library(latex2exp)
  L <- ncol(ppp)
  Methodlist <- c("ipw","reg","dr")
  stat.model <- list(c("T","T","T"),c("M","F","T"),c("T","T","F"),c("M","F","F"))
#############################data pre-processing################################
  rslt.p.value.make <- function(ppp){
    L <- ncol(ppp)
    ppp1 <- rbind(ppp[1,],ppp[1,],ppp[2,],ppp[2,],ppp[3:10,])
    Methodlist <- c("ipw","reg","dr")
    # stat.model <- list(c(T,F),c(F,F),c(F,T))
    stat.model <- list(c("T","T","T"),c("M","F","T"),c("T","T","F"),c("M","F","F"))
    WholeL <- length(stat.model)*length(Methodlist)*L
    rslt.p.value <- data.frame(PValues = rep(0,WholeL),Method = rep('0',WholeL),Model = rep('0',WholeL))
    rslt.p.value$Model <- as.character(rslt.p.value$Model)
    rslt.p.value$Method <- as.character(rslt.p.value$Method)
    t <- 1
    r <- 1
    for (Method in Methodlist) {
      for (k in 1:length(stat.model)){
        model <- switch (paste(stat.model[[k]],collapse = ""),
                         'TTT' = 'True ps; True reg',
                         'MFT' = 'True ps; Wrong reg',
                         'TTF' = 'Wrong ps; True reg',
                         'MFF' = 'Wrong ps; Wrong reg')
        # model <- switch (paste(as.character(stat.model[[k]]),collapse = ""),
        #                  'TRUEFALSE' = 'Strong Null; Unstudentized', 
        #                  'FALSEFALSE' = 'Weak Null; Unstudentized', 
        #                  'FALSETRUE' = 'Weak Null; Studentized') 
        rslt.p.value[t:(t+L-1),1] <- ppp1[r,]
        rslt.p.value[t:(t+L-1),2] <- rep(Method,L)
        rslt.p.value[t:(t+L-1),3] <- rep(model,L)
        t <- t + L
        r <- r + 1
      }
    }
    return(rslt.p.value)
}
#############################data pre-processing################################
  rslt.pdr.value.make <- function(p.values,method){
    if(anyNA(p.values))
      return(p.values)
    else{
      L <- ncol(p.values)
      num <- nrow(p.values)
      return(data.frame(PValues=c(t(p.values[(num-3):num,])),
                          Method=rep(method,4*L),
                          Model=rep(c('True ps; True reg',
                                      'True ps; Wrong reg',
                                      'Wrong ps; True reg',
                                      'Wrong ps; Wrong reg'),each = L)))
    }
  }
#The default mode, plot the p-values for ipw, reg, dr estimator under four scenariors#
if(m=='single'){
  rslt.p.value <- rslt.p.value.make(ppp)
  rslt.p.value$Method <- factor(rslt.p.value$Method,levels = Methodlist)
  plot.p.value <- ggplot(rslt.p.value, aes(x=PValues)) +
    geom_histogram(breaks = 0:20/20,aes(y=..density..),colour="gray77",fill = "gray77") +
    geom_hline(yintercept=1,color="black", linetype="dashed", size = 0.5) +
    theme_set(theme_bw()) + scale_y_continuous(limits = c(0,2),oob = scales::squish) +
    theme(panel.grid.major = element_line(colour=NA),panel.grid.minor = element_blank()) +
    facet_grid(rows = vars(Method), vars(cols = Model), scales = "free_y")+ 
    labs(x = "p-value")
  return(plot.p.value)
}
else if (m=="cmp"){
################compare four testing methods under four scenariors##############
  ppp_dr <- rslt.pdr.value.make(ppp,"PPP (asymp sd)")
  normal_dr <- rslt.pdr.value.make(p.values,"normal (asymp sd)")
  boot_dr <- rslt.pdr.value.make(ppp_boot,"normal (boot sd)")
  bootppp_dr <- rslt.pdr.value.make(bootppp,"PPP (boot sd)")
  rslt_level <- c(ifelse(anyNA(normal_dr),NA,"normal (asymp sd)"),
                  ifelse(anyNA(boot_dr),NA,"normal (boot sd)"),
                  ifelse(anyNA(ppp_dr),NA,"PPP (asymp sd)"),
                  ifelse(anyNA(bootppp_dr),NA,"PPP (boot sd)"))
  rslt_level <- rslt_level[!is.na(rslt_level)]
  rslt.p.value2 <- rbind(ppp_dr,normal_dr,boot_dr,bootppp_dr)
  rslt.p.value2 <- na.omit(rslt.p.value2)
  rslt.p.value2$Method <- factor(rslt.p.value2$Method,levels = rslt_level)

  plot.p.value <- ggplot(rslt.p.value2, aes(x=PValues)) +
    geom_histogram(breaks = 0:20/20,aes(y=..density..),colour="gray77",fill = "gray77") +
    geom_hline(yintercept=1,color="black", linetype="dashed", size = 0.5) +
    theme_set(theme_bw()) + scale_y_continuous(limits = c(0,2),oob = scales::squish) +
    theme(panel.grid.major = element_line(colour=NA),panel.grid.minor = element_blank()) +
    facet_grid(rows = vars(Method), vars(cols = Model), scales = "free_y")+ 
    labs(x = "p-value")
  return(plot.p.value)
}
else if (m=="flip"){
################put two plots together, didn't use in this article##############
  rslt.p.value1 <- rslt.p.value.make(ppp) %>% mutate(DGP = "original")
  rslt.p.value2 <- rslt.p.value.make(p.values) %>% mutate(DGP = "flipped")
  rslt.p.value <- rbind(rslt.p.value1,rslt.p.value2) 
  plot.p.value <- ggplot(rslt.p.value, aes(x=PValues,fill = DGP)) +
    geom_histogram(position="identity" ,breaks = 0:20/20,alpha = 0.8, aes(y=..density..)) +
    scale_fill_manual(values = c("grey22", "grey77")) +
    geom_hline(yintercept=1,color="black", linetype="dashed", size = 0.5) +
    theme_set(theme_bw()) + scale_y_continuous(limits = c(0,2),oob = scales::squish) +
    theme(panel.grid.major = element_line(colour=NA),panel.grid.minor = element_blank()) +
    facet_grid(rows = vars(Method), vars(cols = Model), scales = "free_y")+ 
    labs(x = "p-value")
  plot.p.value
}
else if (m=="real"){
######the posterior predictive value of three kinds of test statistics##########
############################didn't use in this article##########################
  rslt <- data.frame(value = c(ppp,p.values,ppp_boot),stat = rep(rep(c("ipw","reg","dr"),each = params$m),3),class = rep(c("no match","1-1 match","1-2 match"),each = 3*ncol(ppp)))
  rslt$stat <- factor(rslt$stat, levels = c("ipw","reg","dr"))
  rslt$class <- factor(rslt$class, levels = c("no match","1-1 match","1-2 match"))
  ggplot(rslt,aes(x = value)) + scale_x_continuous(limits = c(-3,3),oob = scales::censor) +
  geom_histogram(binwidth = 0.15,aes(y=..density..),colour="gray77",fill = "gray77") +
  facet_grid(class~stat, scales = "fixed") + 
  geom_line(aes(y = ..density..,colour = 'Empirical'),stat = 'density') +
  stat_function(fun = dnorm, args = list(0,1), aes(colour = "normal")) +
  labs(x = TeX("$\\hat{\\tau}^r/se^r$")) + scale_colour_hue("density lines")
}
}

